Dynamic Programming Optimization

Optimal load shedding for Aggregate Induction Motor Load with UPS

Harry G. Kwatny
Edoe Mensah

Department of Mechanical Engineering & Mechanics
Drexel University
Philadelphia PA 19104

1. Introduction

Often, the role of logic in mathematical programming goes unrecognized. Three decades ago,  it was shown that  posing logical inference as a mathematical programming problem opens up the possibility of transfering the computational techniques of mathematical programming to the symbolic world of logic.
For novice decision makers the classical approach of integer programming modeling based on a quantitative paradigm, presents some limitations in introducing qualitative features of modern day complex decision making. In addition to the inability of a non-specialist to interprete purely quantitative results, and its limitations for specification of  real world problems, researchers have developed and continue to develop new algorithms that will enable the modeling of qualitative information often encountered in complex decision support procedures. 
The most significant contribution came from Robert Jeroslow and Paul Williams in their works on the connections between logic and mathematical programming [4]. Granot and Hammer in 1970, Hammer and Peled in 1972 had exploited the analogy between Boolean algebra and 0-1 programming. And later F. Glover in 1975 represents coupling between logical and continuous variables as linear inequalities on binary and continuous variable. Based on these preexisting logic ideas Morari and Bemporad recently devised the Mixed logical Dynamical systems framework for hybrid systems [9]
Qiang Li, Yike Guo and Tetsuo Ida had proposed a modeling language for IP based on the propositional logic and a systematic algorithm for transforming logic formulas into IP formulas using mathematica.
Our goal is to improved the proposed  modeling language and use it for modeling hybrid systems.

Dynamic programming

2. Propositional Logic

Atomic formula: are proposition without logical connectives.
It will state either  an object is a member of a set, or that it is equal to another set.
Valid atomic proposition: 
Tomorrow is tuesday: true if today is monday False if not
2+2=5:False

The syntax is the set of rules defining how atomic formula may be combined to well-formed formulas.
Every atomic formula is a well-formed formula
An alphabet of symbols:(∧,∨,¬,⇒,⇔(,), a,b,..z, a1,b1,..z1,....)
Rules:
Negation: ¬p
Conjunction: p∧q
Disjunction:  p∨q
Implication:  p⇒q
Equivalence:p⇔q
If p and q are well-formed formulas then, so are:(¬p), (¬q), ( p∧q), ( p∨q), (p⇒q), (p⇔q)

It characterize the meaning of the logical proposition with respect to a specified frame. The frame is a set of states through which a system might evolve ( ie hybrid automaton)

A clause is a logical formula with literals in disjunctive form.Ex: x1 ∨ ¬x2∨x3
A formula is in clausal form(CF) or Conjunctive Normal From(CNF) is it a conjunction of clauses. 
Ex:(x1 ∨ ¬x2∨x3)∧(x1 ∨ ¬x2∨¬x3)
The dual of the conjunctive normal form is the disjunctive normal form. it is the disjunction of conjounctive literals: 
Ex: (x1∧ ¬x2)∨(x2∧¬x3)
 A boolean function f(x) is defined by the list of assignments that make it false.This function can be express in CNF by writing a formula that rules out all such assignment.
 Let write a boolean function which is  false when (x1,x2)=(1,1) or (0,0) and true otherwise:(¬x1∨¬x2)∧(x1∨x2) 
 Constraint satisfaction:(satisfiability)
 Any constraint whose variables are two-valued can be interpreted as defining a boolean function f(x) that is false when x is violate the constraint.
 Ex : the 0-1 inequality 3x1+2x2+x3≥4 defines a set a boolean function f(x) that is false when x1=0 or (x2,x3)=(0,0)  and true otherwise.
 The inequality is therefore equivalent to the clause set: { x1 , x2∨x3}
 A dual principle operates for DNF formulas. A DNF equivalent of a boolean function f(x1,x2) defined by (¬x1∨¬x2)∧(x1∨x2), which is true for (x1,x2)=(0,1) and (x1,x2)=(1,0) , can be written as the DNF  formula.  (¬x1∧ x2)∨(x1∧¬x2)
 Neither a CNF nor DNF representation of a boolean function is unique .
 For example the following  CNF formulas represent the same boolean function: (x1∨x2)∧(¬x1∨x2) , (x1∨x2)∧(¬x1∨x2)∧(x2∨x3).
 In general compound proposition will have many different representation. In order to give comparable representations it is convenient to use one of the two common standard representations. CNF, DNF.
 Note: all unsatisfiable CNF formulas represent the same function, namely the constant function f(x)=0 and the DNF formulas represent f(x)=1 
 IP Formulas:

To each propositional variable p, we can associate a decision variable  p , a variable that can assume only two integer values {1,0}. 
(Subscript[δ, p]) = 1 means p=True
(Subscript[δ, p])= 0 means p=False

IP formulas are Integer programming
inequality constraints involving decision variables.
Or:        p ∨ q  ⟷ (Subscript[δ, p]) + (Subscript[δ, q])  ≥1
And:       p ∧ q  ⟷ (Subscript[δ, p]) ≥1,  (Subscript[δ, q])  ≥1
If then:    p ⇒q ⟷ ¬p ∨ q  ⟷1- (Subscript[δ, p]) + (Subscript[δ, q])≥1⟷  (Subscript[δ, q]) ≥ (Subscript[δ, p])
IFF:       p ⇔ q ⟷ (p ⇒ q) ∧ (q  ⇒ p)⟷ (Subscript[δ, q])  ≥ (Subscript[δ, p]) , (Subscript[δ, p])  ≥ (Subscript[δ, q])
XOr:      p ⊕ q ⟷ (p ∨ q )∧(¬ p∨¬ q ) ⟷ (Subscript[δ, p]) + (Subscript[δ, q])  ≥1, (Subscript[δ, p])+(Subscript[δ, q])≤1

3.  IP Functions

Function  ΓFormQ[x]  checks whether an expression is  Γ-form .

ΓFormQ[Γ[_,_]] := True;
ΓFormQ[_] := False;

Function VarQ[x] checks whether x is a propositional variable or an integer variable.

VarQ[x_Symbol]/;x=!=True  ∧x=!=False :=True; 
 VarQ[_]:=False;

Function NegPVarQ[x] checks whether x is a negated propositional variable.

NegPVarQ[¬x_Symbol]/;x=!=True  ∧x=!=False:= True;
NegPVarQ[_]:=False;

Function PLiteralQ[p] checks whether p is a propositional variable or a negated  propositional variable

PLiteralQ[p_]:=VarQ[p]∨NegPVarQ[p];
NotPLiteralQ[p_]:=¬PLiteralQ[p]

Function BoundedQ[x] checks whether a variable x is properly bounded.

BoundedQ[m_Integer≤x_?VarQ≤n_Integer]:=True;
BoundedQ[m_Real≤x_?VarQ≤n_Integer]:=True;
BoundedQ[m_Integer≤x_?VarQ≤n_Real]:=True;
BoundedQ[m_Real≤x_?VarQ≤n_Real]:=True;
BoundedQ[m_Integer<x_?VarQ<n_Integer]:=True;
BoundedQ[m_Real<x_?VarQ<n_Integer]:=True;
BoundedQ[m_Integer<x_?VarQ<n_Real]:=True;
BoundedQ[m_Real<x_?VarQ<n_Real]:=True;
BoundedQ[_]:= False;

Function ConstraintQ[c] checks whether an expression c  is a  constraint . Note that we use ==  to  denote the equality of symbol of (SuperPlus[L]).

ConstraintQ[x_≥y_|x_>y_|x_≤y_ | x_<y_  | x_==y_ ]:= True;
ConstraintQ[_] := False;

Function  PredOrder[x,y] defines a predicate order in a Γ-formula.

PredOrder[x_?VarQ,r_[_,_]] := True;
PredOrder[Not[_],r_[_,_]] := True;
PredOrder[x_,y_] :=False

Function PLiteral2IPvar[e] translates a propositional variable e= p and a formula e=¬p into IP-formulas.

PLiteral2IPvar[x_?VarQ]:=Bounds[[Position[Bounds,x][[1,1]],2]] ;
PLiteral2IPvar[¬x_?VarQ]:=1-Bounds[[Position[Bounds,x][[1,1]],2]] ;

Function Raux[x,y,z] is an auxiliary function of RuleR[x].

Raux[{x___?PLiteralQ,y___ ?NotPLiteralQ},m_,n_]  :=Module[{newvarlist=Map[Unique["d"]&,{y}]},Bounds=Bounds⋃Map[0≤#≤1&,newvarlist];
		Prepend[Map[RuleR,MapThread[Function[{d,a},d>0 \[Implies]a],{newvarlist,{y}}]],(Apply[Plus,newvarlist⋃Map[PLiteral2IPvar,{x}]]≥m  n)]];

Function RuleE[x] defines rule E.

RuleE1[x_]:=x//.{p_⇒q_->¬p∨q };
RuleE2[x_] :=x//.{p_\[Implies]q_->¬p∨q,p_⊕q_->(¬q∨¬p)∧(p∨q),p_⇔q_->( p\[Implies]q)∧(q\[Implies]p)};
RuleE3[x_] := x//.{p_⊕q_->(¬q∨¬p)∧(p∨q)};
RuleE4[x_] := x//.{p3_⇔(p1_⊕p2_)->(p2∨¬p1∨p3)∧(p1∨¬p2∨p3)∧(¬p3∨¬p1∨¬p2)∧(¬p3∨p1∨p2),
			 p1_⊕p2_⇒p3_->(p2∨¬p1∨p3)∧(p1∨¬p2∨p3),
			 p3_⇒p1_⊕p2_->(¬p3∨¬p1∨¬p2)∧(¬p3∨p1∨p2),     p_⊕q_->(¬q∨¬p)∧(p∨q)};

Function  R[x] defines rule for Conjunctive Normal Form.

R1[x_] := x//.{p_⊕q_->(¬p∨¬q)∧(p∨q)};
R2[x_]:=x//.{A_⇒B_->(¬A∨B),A_⇔B_->(¬A∨B)∧(A∨¬B),(A_∨¬A_)->True,(¬A_∨A_)->True};
R3[x_]:=x//.{¬(A_∧B_)->(¬A∨¬B),¬(A_∨B_)->(¬A∧¬B)};
R4[x_]:=x//.{(A_∧B_)∨C_->(A∨C)∧(B∨C),C_∨(A_∧B_)->(C∨A)∧(C∨B)};
R5[x_]:=x//.{(A_∨¬A_)->True,(¬A_∨A_)->True};
R6[x_]:=x//.{(A_∨A_)->A};
CNF[Predicate_∧x_]:=Predicate∧CNF[R2[R1[RuleE1[RuleE2[x]]]]]
CNF[f_]:=Not[LogicalExpand[Not[R2[R1[f]]]]]//.{Not[a_Or]:>And@@(Not/@List@@a),Not[a_And]:>Or@@(Not/@List@@a)}//.A_⇒B_->¬A∨B

Function  E[x] defines rule for Disjunctive Normal Form.

E1[x_] := x//.{p_⊕q_->(¬p∨¬q)∧(p∨q)};
E2[x_]:=x//.{(A_⇒B_)->(¬A∨B),((A_)⇒B_)->(¬A∨B),((A_)⇒B_)->¬A∨B,A_⇒B_->¬A∨B,A_⇔B_->(¬A∨B)∧(A∨¬B)};
E3[x_]:=x//.{¬(A_∧B_)->(¬A∨¬B),¬(A_∨B_)->(¬A∧¬B)}
E4[x_]:=x//.{(A_∨B_)∧C_->(A∧C)∨(B∧C),C_∧(A_∨B_)->(C∧A)∨(C∧B)};
E5[x_]:=x//.{(A_∨¬A_)->True,(¬A_∨A_)->True};
E6[x_]:=x//.{(A_∨A_)->A};
DNF[spec_]:=Composition[E6,E5,E4,E3,E2,E1][spec]

Function  RuleT[x] defines rule T.

Clear[RuleT]

RuleT[x1_∨x2_]:=Γ[1,{RuleT[x1],RuleT[x2]}];
RuleT[x1_∧x2_]:=Γ[2,{RuleT[x1],RuleT[x2]}];
(*
RuleT[x1_∧x2_∧x3_]:=Γ[3,{RuleT[x1],RuleT[x2],RuleT[x3]}];RuleT[x1_∧x2_∧x3_∧x4_]:=Γ[4,{RuleT[x1],RuleT[x2],RuleT[x3],RuleT[x4]}];
RuleT[x1_∧x2_∧x3_∧x4_∧x5_]:=Γ[5,{RuleT[x1],RuleT[x2],RuleT[x3],RuleT[x4],RuleT[x5]}];
RuleT[x1_∧x2_∧x3_∧x4_∧x5_∧x6_]:=Γ[6,{RuleT[x1],RuleT[x2],RuleT[x3],RuleT[x4],RuleT[x5],RuleT[x6]}];
RuleT[x1_∧x2_∧x3_∧x4_∧x5_∧x6_∧x7_]:=Γ[7,{RuleT[x1],RuleT[x2],RuleT[x3],RuleT[x4],RuleT[x5],RuleT[x6],RuleT[x7]}];
RuleT[x1_∧x2_∧x3_∧x4_∧x5_∧x6_∧x7_∧x8_]:=Γ[8,{RuleT[x1],RuleT[x2],RuleT[x3],RuleT[x4],RuleT[x5],RuleT[x6],RuleT[x7],RuleT[x8]}];
RuleT[x1_∧x2_∧x3_∧x4_∧x5_∧x6_∧x7_∧x8_∧x9_]:=Γ[9,{RuleT[x1],RuleT[x2],RuleT[x3],RuleT[x4],RuleT[x5],RuleT[x6],RuleT[x7],RuleT[x8],RuleT[x9]}];
RuleT[x1_∧x2_∧x3_∧x4_∧x5_∧x6_∧x7_∧x8_∧x9_∧x10_]:=Γ[10,{RuleT[x1],RuleT[x2],RuleT[x3],RuleT[x4],RuleT[x5],RuleT[x6],RuleT[x7],RuleT[x8],RuleT[x9],RuleT[x10]}];
RuleT[x1_∧x2_∧x3_∧x4_∧x5_∧x6_∧x7_∧x8_∧x9_∧x10_∧x11_]:=Γ[11,{RuleT[x1],RuleT[x2],RuleT[x3],RuleT[x4],RuleT[x5],RuleT[x6],RuleT[x7],RuleT[x8],RuleT[x9],RuleT[x10],RuleT[x11]}]; 
*)
RuleT[¬x_]:=¬RuleT[x] ;
RuleT[x_≠y_]:=¬RuleT[x==y];
RuleT[atleast[m_,S_]]:=Γ[m,Map[RuleT,S]];
RuleT[atmost[m_,S_]]:=Γ[Length[S]-m,Map[RuleT,Map[Not,S]]];
RuleT[none[S_]]:=Γ[Length[S],Map[RuleT,Map[Not,S]]];
RuleT[x_?VarQ]:=PVarDecl[x];
RuleT[x_?ConstraintQ]:=IPVarDecl[x];
RuleT[True]:=Throw[{True,Bounds}];
RuleT[False]:=Throw["Wrong specification"];
RuleT[exactly[m_,S_]]:=Γ[2,{RuleT[atleast[m,S]],RuleT[atmost[m,S]]}];

Function RuleN [x] defines rule N.

RuleN[Γ[m_,S_]]:=Γ[m,Map[RuleN,S]];
RuleN[¬Γ[m_,S_]]:=Γ[Length[S]-m+1,Map[RuleN,Map[Not,S]]];
RuleN[¬(x1_==x2_)]:=Γ[1,{x1>x2,x1<x2}];
RuleN[x_]:=x;

Function  RuleF[x] defines rule  F.

RuleF[x_]:=x //.{RF1,RF2}
RF1= Γ [1,{x1___,Γ[1,{x2___}],x3___}] :> Γ[1,{x1,x2,x3}];
RF2 = Γ[m_,{x1___,Γ[n_,S_],x2___}]/;m==Length[{x1,x2}]+1∧n==Length[S] :>  Γ[n+m-1,S⋃{x1,x2}];

Function  RuleR[x] defines rule R.

RuleR[r_[x1_,x2_]]:=r[x1,x2];
RuleR[p_?PLiteralQ]:=PLiteral2IPvar[p]⩾1;
RuleR[Γ[m_,S_]]/;m==Length[S]:=Map[RuleR,S];
RuleR[Γ[m_, S_]] /;m<Length[S]:= Raux[Sort[S,PredOrder],m,1];
RuleR[d_>0\[Implies]Γ[m_,S_]]/;m==Length[S]:=Map[RuleR,Map[(d>0\[Implies]#)&,S]];
RuleR[d_>0\[Implies]Γ[m_,S_]]/;m<Length[S]:=Raux[Sort[S,PredOrder],m,d];
RuleR[d_>0\[Implies]p_?PLiteralQ]:=PLiteral2IPvar[p]⩾d;
RuleR[d_>0\[Implies]x1_⩾x2_]:= Module[{L},(Subscript[(x1-x2 ⩾ L(1-d )), L])];
RuleR[d_>0\[Implies]x1_>x2_]:=Module[{L},(Subscript[(x1-x2 > (L-1) (1-d )), L])];RuleR[d_>0\[Implies]x1_⩽x2_]:=Module[{U},(((x1-x2 ⩽U (1-d )))^U)];
RuleR[d_>0\[Implies]x1_<x2_]:=Module[{U},(((x1-x2 <(U+1) (1-d )))^U)];
RuleR[d_>0\[Implies]x1_==x2_]:={RuleR[d>0\[Implies]x1⩾x2],RuleR[d>0\[Implies]x1⩽x2]};

Function Lowerbound[A] computes the  lowerbound of  constraint A with  respects to the  given   Bounds, where  Bounds is specifed as  a list of terms of the form a≤x≤b.   
Lowerbound uses  auxiliary  function Lb.

Lowerbound[A_?AtomQ ]:=Lb[A ];
Lowerbound[A_ ]:= Map[Lb,A];
Lb[term:a_*  x_ ] /; a>0 :=Module[{i=Position[Bounds,_≤x≤_][[1,1]]},term/.{Bounds[[i,2]]->Bounds[[i,1]]}];
Lb[term:a_ * x_ ] /; a<0  :=Module[{i=Position[Bounds,_≤x≤_][[1,1]]},term/.{Bounds[[i,2]]->Bounds[[i,3]]}];
Lb[x_?VarQ ]  :=Bounds[[Position[Bounds,_≤x≤_][[1,1]],1]]
Lb[a_ ] :=a

Function Upperbound[A] computes the  lowerbound of  constraint A with  respects to the  given   Bounds, where  Bounds is specifed as  a list of terms of the form a≤x≤b.   
Upperbound uses  auxiliary  function Lb.

Upperbound[A_]:= Map[Ub,A]; 
Upperbound[A_?AtomQ]:=Ub[A ];
Ub[term:a_*x_ ] /; a>0 :=Module[{i=Position[Bounds,_≤x≤_][[1,1]]},term/.{Bounds[[i,2]]->Bounds[[i,3]]}];
Ub[term:a_* x_ ] /; a<0  :=Module[{i=Position[Bounds,_≤x≤_][[1,1]]},term/.{Bounds[[i,2]]->Bounds[[i,1]]}];
Ub[x_?VarQ ] :=Bounds[[Position[Bounds,_≤x≤_][[1,1]],3]]; 
Ub[a_ ] := a;

Function RRC[x] (Remove Redundant Constraints) is defined below using the following    auxiliary functions:  
RCaux  checks whether a given constraint is redundant and outputs {} if redundant or else  a simplified constraint after substituting the    value of the Bounds.

RRC[A_]:= Map[RCaux,A];
RCaux[ (((A_≤B_))^(U_)) ]:=Module[{Uw=Upperbound[A]},
		If[Uw>0,A≤(B/.{U->Uw}),{}]];
RCaux[(((A_<B_) )^(U_))]:=
	Module[{Uw=Upperbound[A]},
		If[Uw≥0,A<(B/.{U->Uw}),{}]];
RCaux[ (Subscript[(A_≥B_) , L_])] :=Module[{Lw=Lowerbound[A]},
		If[Lw<0,A≥(B /.{L->Lw}),{}]];
RCaux[(Subscript[(A_>B_) , L_])]:=
	Module[{Lw=Lowerbound[A]},
		If[Lw≤0,A>(B/.{L->Lw}),{}]];
RCaux[A_]:=A;

Function MakeBounds[x] checks whether the input bound of each integer variable is right and   generates a decision variable and its bound for each  propositional    variable.

MakeBounds[ x_?BoundedQ ]:= x;
MakeBounds[x_?VarQ]:=0≤(Subscript[Pvar, x])≤1;

Founction PVarDecl[e] checks  whether each propositional variable is properly declared.

PVarDecl[v_]:=If[MemberQ[Bounds,v,{2,3}],v,Print["Propositional Variable ",  v," not declared."]; Throw[Bounds]];

Founction IPVarDecl[e] checks  whether each integer variable  is properly declared.

IPVarDecl[v_?VarQ]:=If[MemberQ[Bounds,v,2],v,Print["Integer variable ",  v," not declared."]; 
		Throw[Bounds]];
IPVarDecl[t:f_[__]]:=Map[IPVarDecl,t];
IPVarDecl[x_Integer]:=x;
IPVarDecl[x_Real]:=x;
IPVarDecl[x_]:= Module[{},Print["Constraint ",x," is not well formed"];Throw[Bounds]]

Founction IneqSimplify[e] simplifies expressions.

IneqSimplify[a_> b_] := Simplify[a -b] >0;
IneqSimplify[a_≥ b_] := Simplify[a -b]  ≥0;
IneqSimplify[a_< b_] := Simplify[a -b] <0;
IneqSimplify[a_≤ b_] := Simplify[a -b]  ≤0;

Founction Options[GenIP]  specifies that the generated IP-formulas will be simplified.

Options[GenIP]={IPFormSimplify->True};

Main function GenIP.  GenIP needs three arguments spec , decl and opts.  
spec is  a logical specification, 
decl is a list of propositional variables and  the bounds of  integer variables of the form  a ≤x≤b, and 
opts is an option to specify whether the generated IP-formulas will be simplified or not with IPFormSimplify->True  or IPFormSimplify->False.

GenIP[spec_,decl_,opts___]:= Catch[Block[{Pvar=δ ,Bounds=Map[MakeBounds,decl]} ,	 
  With[{ipform=Flatten[{Composition[RRC,Flatten,RuleR,RuleF,RuleN,RuleT,RuleE2,RuleE1,RuleE4,RuleE3][spec]}]},  If[IPFormSimplify/.{opts}/.Options[GenIP],Map[IneqSimplify,ipform],ipform,ipform]⋃Bounds]]]

Specification of nice

Let p be a boolean expression, then nice[p] gives p in "nice" form (i.e. "||" is replaced by "+", "&&" is replaced by "·" and "!" is replaced by "(OverBar[□])"):

nice[p_]:=p/.{And->Times,Or->Plus,Not->OverBar};

nice[(a||b)&&(c||!d)||!a]

unnice[p_]:=p/.{Times->And,Plus->Or,OverBar->Not};

4. Dynamic Programming Functions

The main function is OptimalPolicy which calls the function Findmin. Given the current state, FindMin computes the optimal control that will be applied over the next time interval.

FindMin

Define binary and real variables: BinaryVars, RealVars
Split Inequalities into eqBinary, eqReal, binary equations contain only binary variables, real equations can contain both binary and real variables
GIVEN RealStateVars, BinaryStateVars, BinaryVars, RealVars, eqBinary, eqReal, J = cost as (interpolation) function of  (RealStateVars, BinaryVars, RealVars, AltBinaryState)
Use Reduce to get all feasible binary combinations. Note: In this case there are 2^8=256 possible combinations of binary values, but only 8 are feasible
Use Reduce to solve for real variables for each feasible combination of binary variables - the result will be unique for each feasible combination of binary variables
Evaluate J for each pair of binary & real & select minimum

Clear[FindMin];
FindMin[ListSolsBinary_List,EqReal_List,BinaryVars_List,RealVars_List,RealStateVars_List,BinaryStateVars_List,x_List,AltBinaryState_,NextBinaryState_,J_]:=Module[{StateReplace,SolsBinary,SolsReal,Sols,Cost,Ind,A,RR,SolAlt},
StateReplace=Inner[Rule,RealStateVars,Drop[x,-1],List];
RR=True;
Map[(RR=RR&&#)&,EqReal];
SolsReal=Map[Resolve[Exists[RealVars,Rationalize[(RR//.StateReplace)]//.#],Reals]&,ListSolsBinary[[x[[-1]]]]];
(* SolsBinary is solved and the result is substituted in the EqReal. Similarily EqReal now contained only real variable is ratinalized and reduced *)
Ind=Position[SolsReal,True];
(* Identify positions where there are real solution.  *)
SolsReal=N[Map[Reduce[Rationalize[(EqReal//.StateReplace)]//.#,RealVars,Reals]&,ListSolsBinary[[x[[-1]]]][[Flatten[Ind]]]]];
   A=ListSolsBinary[[x[[-1]]]][[Flatten[Ind]]];
(*  Take the corrsponding binary solutions *)
Sols=MapThread[Join[#1,#2]&,{A,Map[Solve[#,RealVars][[1]]&,SolsReal]}];
 (* Pair up the valid binary and real solutions *)
(* Set up a rule to switch to the alternate binary state representation *)
 SolAlt=Map[{AltBinaryState->Position[NextBinaryState/.#,1][[1]]}&,Sols];
(* SolAlt select the NextBinaryState (if any) for which a transition occurs, ie it binary value is 1  *)
Cost=Flatten[MapThread[(((J//.StateReplace)//.#1)//.#2)&,{Sols,SolAlt}]];
(* A mapping of J (cost) with the Sols, after SolAlt  is substituted  *)
   Ind=Position[Cost,Min[Cost]]; (* assign a position to the minimum cost*)
  (*  If[Length[Ind]>1,Print["Warning: Multiple minima found! Indices = "<>ToString[Ind]]]; *)
   {Cost[[Ind[[1]]]],Sols[[Ind[[1]]]]}
(* Select the cost value and the solution Sols corresponding to the minimum cost*)
]

Clear[RealSolutions]
RealSolutions[ListSolsBinary_List,EqReal_List,BinaryVars_List,RealVars_List,RealStateVars_List,BinaryStateVars_List,x_List]:=Module[{StateReplace,SolsBinary,SolsReal,Sols,Ind,A,RR},
(* Given full state x, and list of binary solutions, find real solutions. RealSolutions will be mapped over S. *)
StateReplace=Inner[Rule,RealStateVars,Drop[x,-1],List]; (* Pull out real state from x *)
RR=True;
Map[(RR=RR&&#)&,EqReal]; (* Convert real equation list to "&&" *)
SolsReal=Map[Resolve[Exists[RealVars,Rationalize[(RR//.StateReplace)]//.#],Reals]&,ListSolsBinary[[x[[-1]]]]];
(* For each binary solution associated with given discrete state, substitute current real state into EqReal and check if next real state exists using Resolve *)
Ind=Position[SolsReal,True];
(* Identify positions where there are real solutions.  *)
SolsReal=N[Map[Reduce[Rationalize[(EqReal//.StateReplace)]//.#,RealVars,Reals]&,ListSolsBinary[[x[[-1]]]][[Flatten[Ind]]]]]; (* use reduce to find real solutions if where they exist *)
   A=ListSolsBinary[[x[[-1]]]][[Flatten[Ind]]];
(*  Take the corrsponding binary solutions *)
Sols=MapThread[Join[#1,#2]&,{A,Map[Solve[#,RealVars][[1]]&,SolsReal]}];
 (* Pair up the valid binary and real solutions *)
Sols
]
Clear[AllSolutions]
AllSolutions[ListSolsBinary_List,EqReal_List,BinaryVars_List,RealVars_List,RealStateVars_List,BinaryStateVars_List,S_List]:=Map[RealSolutions[ListSolsBinary,EqReal,BinaryVars,RealVars,RealStateVars,BinaryStateVars,#]&,S]

Clear[FindMin2];
FindMin2[ListSolsAll_List,x_List,AltBinaryState_,NextBinaryState_,J_]:=Module[{StateReplace,Cost,Ind,SolAlt},
StateReplace=Inner[Rule,RealStateVars,Drop[x,-1],List]; (* Pull out real state from x *)
(* Set up a rule to switch to the alternate binary state representation *)
    SolAlt=Map[{AltBinaryState->Position[NextBinaryState/.#,1][[1]]}&,ListSolsAll];
(* Compute cost associated with each admissible mode change *)
Cost=Flatten[MapThread[(((J//.StateReplace)//.#1)//.#2)&,{ListSolsAll,SolAlt}]];
Ind=Position[Cost,Min[Cost]]; (* identify positions of minimum cost*)
  (*  If[Length[Ind]>1,Print["Warning: Multiple minima found! Indices = "<>ToString[Ind]]]; *)
   {Cost[[Ind[[1]]]],ListSolsAll[[Ind[[1]]]]}
(* Select the cost value and the solution Sols corresponding to the minimum cost*)
]

Auxilliary Functions

The problem real variables are always known from the problem definition. However, the bianry variables need to be identified after the IP fromulas are determined. Here are two functions that generate the binary variables, the binary equations and the real equations and a simple one that generates the real equations after binary equations are identified.

BinaryEquations[IPFormulas_List,RealVars_List]:=Module[{DropReals,AA=IPFormulas},
DropReals[AA_List,i_]:=AA[[Flatten[Position[Map[(Cases[Variables[Level[#,1]],i])&,AA],{}]]]];
Map[(AA=DropReals[AA,#])&,RealVars];
AA
];

RealEquations[IPFormulas_,EqBinary_]:=Complement[IPFormulas,EqBinary];

BinaryVariables[IPFormulas_,RealVars_]:=Module[{AA=IPFormulas},Complement[Union[Flatten[Map[(Variables[Level[#,1]])&,AA]]],RealVars]];

ControlRules[ControlLogic_,StateReplaceList_,DiscreteState_,InputEvents_]:=Module[{ControlIP,F,vars,ControlEqns,CombinedList,ControlReplacements},
ControlIP=Map[GenIP[#,Join[DiscreteState,InputEvents]]&,ControlLogic];
F[x_]:=Union[Flatten[Map[Variables[Level[#,1]]&,x]]];
vars=Map[F[#]&,MapThread[#1/.#2&,{ControlIP,StateReplaceList}]];
ControlEqns=MapThread[Reduce[#1/.#2,#3,Integers]&,{ControlIP,StateReplaceList,vars}];
CombinedList=MapThread[{#1,#2}&,{ControlEqns,vars}];
ControlReplacements=Map[Solve[#[[1]],#[[2]]]&,CombinedList];
ControlReplacements
]

BinarySolutions[EqBinary_,StateReplaceList_,ControlReplacements_]:=Module[{MyEqns,F,Ans,G,Sols,H},
MyEqns=MapThread[(EqBinary//.#1)//.#2&,{StateReplaceList,ControlReplacements}];
F[x_]:=Map[Reduce[#,Union[Flatten[Map[Variables[Level[#,1]]&,#]]],Integers]&,x];
Ans=Map[F[#]&,MyEqns];
(* x corresponds to discrete state results (Ans), y to control (ControlReplacements) *)
G[x_,y_]:=Module[{},
H[k_]:=Module[{DDD,EEE},
DDD=Map[Flatten[Solve[#,Variables[Level[#,2]]]]&,Level[x[[k]],1]];
EEE=Map[Join[#,y[[k]]]&,DDD]
];
Map[H[#]&,Range[Length[y]]]
];
MapThread[Flatten[G[#1,#2],1]&,{Ans,ControlReplacements}]
]

Example (3- bus induction motor)

ControlLogic={q1⇒((¬s1∨s1)∧¬s2),
(q2⇒(¬s2∨s1∨¬(¬s2∨s1))∧atmost[1,{¬s2,s1}]),
q3⇒((s2∨¬s2)∧¬s1)};
eqBinary={-1+d7+d8≥0,1-(Subscript[δ, q1])-(Subscript[δ, q2])-(Subscript[δ, q3])≥0,-1+(Subscript[δ, q1])+(Subscript[δ, q2])+(Subscript[δ, q3])≥0,d4-(Subscript[δ, qq1])≥0,d5-(Subscript[δ, qq2])≥0,d6-(Subscript[δ, qq3])≥0,1-(Subscript[δ, qq1])-(Subscript[δ, qq2])-(Subscript[δ, qq3])≥0,-1+(Subscript[δ, qq1])+(Subscript[δ, qq2])+(Subscript[δ, qq3])≥0,1-(Subscript[δ, q1])+(Subscript[δ, qq1])-(Subscript[δ, s1])≥0,1-(Subscript[δ, q2])+(Subscript[δ, qq1])-(Subscript[δ, s1])≥0,-(Subscript[δ, q1])+(Subscript[δ, qq2])+(Subscript[δ, s1])≥0,1-(Subscript[δ, q3])+(Subscript[δ, qq2])-(Subscript[δ, s2])≥0,1-(Subscript[δ, q2])+(Subscript[δ, qq2])+(Subscript[δ, s1])-(Subscript[δ, s2])≥0,-(Subscript[δ, q2])+(Subscript[δ, qq3])+(Subscript[δ, s2])≥0,-(Subscript[δ, q3])+(Subscript[δ, qq3])+(Subscript[δ, s2])≥0,0≤d4≤1,0≤d5≤1,0≤d6≤1,0≤d7≤1,0≤d8≤1,0≤(Subscript[δ, q1])≤1,0≤(Subscript[δ, q2])≤1,0≤(Subscript[δ, q3])≤1,0≤(Subscript[δ, qq1])≤1,0≤(Subscript[δ, qq2])≤1,0≤(Subscript[δ, qq3])≤1,0≤(Subscript[δ, s1])≤1,0≤(Subscript[δ, s2])≤1};
StateReplaceList={{(Subscript[δ, q1])->1,(Subscript[δ, q2])->0,(Subscript[δ, q3])->0},{(Subscript[δ, q1])->0,(Subscript[δ, q2])->1,(Subscript[δ, q3])->0},{(Subscript[δ, q1])->0,(Subscript[δ, q2])->0,(Subscript[δ, q3])->1}};
ControlReplacements=ControlRules[ControlLogic,StateReplaceList,{q1,q2,q3},{s1,s2}]

BinTest=BinarySolutions[eqBinary,StateReplaceList,ControlReplacements];

Length[Flatten[BinTest]]

AA={V1==0.5009756106837937` (Sqrt[(((4.`+(s^2)) (0.01556420233463035`+0.2490272373540856` s+(s^2)))/(((0.0625`+0.015625` (s^2)))^2))]) V2 (1-η),3-d10-V1>0,1-d10+V1>0,-1+d10+d9≥0,-2 d9+V1≥0,-d10+V2≥0,1-(Subscript[δ, q1])-(Subscript[δ, q2])-(Subscript[δ, q3])≥0,-1+(Subscript[δ, q1])+(Subscript[δ, q2])+(Subscript[δ, q3])≥0,1-(Subscript[δ, qq1])-(Subscript[δ, qq2])-(Subscript[δ, qq3])≥0,-1+(Subscript[δ, qq1])+(Subscript[δ, qq2])+(Subscript[δ, qq3])≥0,1-(Subscript[δ, q1])+(Subscript[δ, qq1])-(Subscript[δ, s1])≥0,1-(Subscript[δ, q2])+(Subscript[δ, qq1])-(Subscript[δ, s1])≥0,-(Subscript[δ, q1])+(Subscript[δ, qq2])+(Subscript[δ, s1])≥0,1-(Subscript[δ, q3])+(Subscript[δ, qq2])-(Subscript[δ, s2])≥0,1-(Subscript[δ, q2])+(Subscript[δ, qq2])+(Subscript[δ, s1])-(Subscript[δ, s2])≥0,-(Subscript[δ, q2])+(Subscript[δ, qq3])+(Subscript[δ, s2])≥0,-(Subscript[δ, q3])+(Subscript[δ, qq3])+(Subscript[δ, s2])≥0,-2+d10+V2≤0,0≤d10≤1,0≤d9≤1,0≤s≤1,0≤V1≤2,0≤V2≤2,0≤η≤1,0≤(Subscript[δ, q1])≤1,0≤(Subscript[δ, q2])≤1,0≤(Subscript[δ, q3])≤1,0≤(Subscript[δ, qq1])≤1,0≤(Subscript[δ, qq2])≤1,0≤(Subscript[δ, qq3])≤1,0≤(Subscript[δ, s1])≤1,0≤(Subscript[δ, s2])≤1}

eqBinary=BinaryEquations[AA,{V1,V2,s,SlipEq,η}]

eqReal={V1==0.5009756106837937` (Sqrt[(((4.`+(s^2)) (0.01556420233463035`+0.2490272373540856` s+(s^2)))/(((0.0625`+0.015625` (s^2)))^2))]) V2 (1-0.4` (Subscript[δ, qq2])-0.8` (Subscript[δ, qq3])),3-d10-V1>0,1-d10+V1>0,-2 d9+V1≥0,-d10+V2≥0,-2+d10+V2≤0,0≤s≤1,0≤V1≤2,0≤V2≤2,0≤0.4` (Subscript[δ, qq2])+0.8` (Subscript[δ, qq3])≤1};
varReal={V1,V2,s,SlipEq,η};
varBinary=BinaryVariables[AA,varReal]

BinTest=BinarySolutions[eqBinary,StateReplaceList,ControlReplacements]

RealSolutions[BinTest,eqReal,varBinary,varReal,{s},{(Subscript[δ, q1]),(Subscript[δ, q2]),(Subscript[δ, q3])},{0.5,2}]

ContinuousS=Flatten[Outer[List,Range[0.1,0.5,0.1]],1];
DiscreteS=Range[3];
S=Map[Flatten[#]&,Flatten[Outer[List,ContinuousS,DiscreteS,1],1]]

Length[S]

ASols=AllSolutions[BinTest,eqReal,varBinary,varReal,{s},{(Subscript[δ, q1]),(Subscript[δ, q2]),(Subscript[δ, q3])},S]

Length[ASols]

Length[Flatten[ASols]]

Length[Flatten[BinTest]]

BinTest[[1]]

ASols[[1]]

BinTest[[2]]

ASols[[2]]

BinTest[[3]]

ASols[[3]]

OptimalPolicy

For now, we assume that we seek binary-valued controls, the optimal policy is the control vector, ustar, determined as a tabular function over the state space, S.
The state is composed of the real state variables, RealStateVars, and the discrete state variables. RealStateVars is simply a vector of length n, where n is the dimension of the real state space, X. We use two different representations of the discrete state space. Suppose there are m modes of the system, labeled (Subscript[q, 1]),...,(Subscript[q, m]). One representation of the discrete state is a list (vector) δ={(Subscript[δ, (Subscript[q, 1])]),..., (Subscript[δ, (Subscript[q, m])])} where each (Subscript[δ, (Subscript[q, i])])is either 0 if the system is not in that state or 1 if it is. Thus there is precisly one element equal to one and all of the others are zero. We denote the space of all δ's, Δ. The state space is Σ=X×Δ. 
The other representation is simply {i} where i is the index such that (Subscript[δ, (Subscript[q, i])])=1. Let Z denote the integers 1,2,... The state space S is the product space S=X×Z. It is easy to see that S⇔Σ. Clearly, given δ, i=Position[δ,1]. or given i, construct δ={(Subscript[δ, (Subscript[q, 1])])=0,...,(Subscript[δ, (Subscript[q, i])])=1,.., (Subscript[δ, (Subscript[q, m])])=0}.
To use OptimalPolicy it is necessary to define the following arguments:
     eqBinary, eqReal, BinaryVars, RealVars, BinaryStateVars, NextBinaryState, RealStateVars, DynamicsRHS, BinaryControlVars, S, GN, G, N

ContinuousS=Outer[List,Range[0.1,1.9,0.01]];
DiscreteS=IdentityMatrix[2];
DiscreteS=Range[2];
State space,S=X×Z = Map[Flatten[#]&,Flatten[Outer[List,ContinuousS,DiscreteS,1],1]];
Integrated cost, G[BinaryVars_List,RealVars_List,StateVars_List,n_Integer]:=(RealVars[[1]]/(1/n))^10/m;
Terminal cost, GN[StateVars_List]:=(0.5 (StateVars[[1]]-qbar)^2);
BinaryStateVars={(Subscript[δ, q1]),(Subscript[δ, q2])};
NextBinaryState={(Subscript[δ, qq1]),(Subscript[δ, qq2])};
BinaryControlVars={(Subscript[δ, s])};
RealStateVars={q};
RealVars={i,z};
AllRealVars={i,q,z};
BinaryVars=Complement[BinaryVariables[AA,AllRealVars],BinaryStateVars];
DynamicsRHS={z};
Number of steps, N

Control=OptimalPolicy[eqBinary,eqReal,BinaryVars,RealVars,BinaryStateVars,NextBinaryState,RealStateVars,DynamicsRHS,BinaryControlVars,S,GN,G,N];

Clear[OptimalPolicy];
OptimalPolicy[EqBinary_List,EqReal_List,ControlLogic_,BinaryVars_List,RealVars_List,BinaryStateVars_List,BooleanStateVars_List,NextBinaryState_List,RealStateVars_List,DynamicsRHS_,BinaryControlVars_List,BooleanControlVars_List,S_,GN_,G_,n_Integer]:=Module[{SR,Data,Jstar,StateVars,NextState,AltBinaryState,sol1,OptSol,Delta,del,StateReplaceList,ControlReplacements,ListSolsBinary,ustar,ustarList={}},
(* Construct data table in the form {{x1,x2,..,xn,Jstar} .....} *)
Data=Map[Join[S[[#]],{Map[GN[#]&,S][[#]]}]&,Range[Length[S]]];
(* Data = a list of the states S paired with corresponding values of GN (the terminal cost) *)
(* Create function Jstar(x1,x2,...,xn) *)
Jstar=Interpolation[Data,InterpolationOrder -> 1];
(* BINARY Solutions - solve binary IP's for all discrete states *)
(* Create an array, Delta, of elements del[i] of the length of BinaryStateVars*)
 Delta=Array[del,Length[BinaryStateVars]];
(* initialization of all elements del to 0 *)
Map[(del[#]=0)&,Range[Length[Delta]]];
StateReplaceList=Map[Inner[Rule,BinaryStateVars,ReplacePart[Delta,1,#],List]&,Range[Length[BinaryStateVars]]];
ControlReplacements=ControlRules[ControlLogic,StateReplaceList,BooleanStateVars,BooleanControlVars];
ListSolsBinary=BinarySolutions[EqBinary,StateReplaceList,ControlReplacements];
(* MAIN Loop, Iterate Backward in Time *)  
For[r=n,r≥1,r--,
   StateVars=Join[RealStateVars,BinaryStateVars]; (* The set of all state variables *)
   NextState=Join[DynamicsRHS,{AltBinaryState}];
(* The set of dynamic variable and NextBinaryState (if any) for which a transition occurs, ie it binary value is 1*)
   sol1[x_ ]:=FindMin[ListSolsBinary,EqReal,BinaryVars,RealVars,RealStateVars,BinaryStateVars,x,AltBinaryState,NextBinaryState, G[BinaryVars,RealVars,StateVars,r]+Apply[Jstar,NextState]];
(* the function sol1 call the function FindMin given it full state x *)
   OptSol=Map[sol1[#]&,S];(*  map sol1 - FindMin- over S( continuous and discrete) and find the minimum  at each state *)
Print["Finished FindMin. r = "];Print[r];
   Data=Map[Join[S[[#]],{OptSol[[Range[Length[S]],1]][[#]][[1]]}]&,Range[Length[S]]];
   Jstar=Interpolation[Data,InterpolationOrder -> 1];
   ustar=Map[Inner[Rule,BinaryControlVars,Flatten[BinaryControlVars//.OptSol[[#,2]]],List]&,Range[Length[S]]];
Print[ustar];
ustarList=Join[ustarList,{ustar}];
];
Simplify[ustarList]
]

Clear[OptimalPolicy2];
OptimalPolicy2[EqBinary_List,EqReal_List,ControlLogic_,BinaryVars_List,RealVars_List,BinaryStateVars_List,BooleanStateVars_List,NextBinaryState_List,RealStateVars_List,DynamicsRHS_,BinaryControlVars_List,BooleanControlVars_List,S_,GN_,G_,n_Integer]:=Module[{SR,Data,Jstar,StateVars,NextState,AltBinaryState,sol1,OptSol,Delta,del,StateReplaceList,ControlReplacements,ListSolsBinary,ustar,ustarList={}},
(* Construct data table in the form {{x1,x2,..,xn,Jstar} .....} *)
Data=Map[Join[S[[#]],{Map[GN[#]&,S][[#]]}]&,Range[Length[S]]];
(* Data = a list of the states S paired with corresponding values of GN (the terminal cost) *)
(* Create function Jstar(x1,x2,...,xn) *)
Jstar=Interpolation[Data,InterpolationOrder -> 1];
(* BINARY Solutions - solve binary IP's for all discrete states *)
(* Create an array, Delta, of elements del[i] of the length of BinaryStateVars*)
 Delta=Array[del,Length[BinaryStateVars]];
(* initialization of all elements del to 0 *)
Map[(del[#]=0)&,Range[Length[Delta]]];
StateReplaceList=Map[Inner[Rule,BinaryStateVars,ReplacePart[Delta,1,#],List]&,Range[Length[BinaryStateVars]]];
ControlReplacements=ControlRules[ControlLogic,StateReplaceList,BooleanStateVars,BooleanControlVars];
ListSolsBinary=BinarySolutions[EqBinary,StateReplaceList,ControlReplacements];
ListSolsAll=AllSolutions[ListSolsBinary,EqReal,BinaryVars,RealVars,RealStateVars,BinaryStateVars,S];
Print["Toatal number of solutions = "<>ToString[Length[Flatten[ListSolsAll]]]];
(* MAIN Loop, Iterate Backward in Time *)  
For[r=n,r≥1,r--,
   StateVars=Join[RealStateVars,BinaryStateVars]; (* The set of all state variables *)
   NextState=Join[DynamicsRHS,{AltBinaryState}];
(* The set of dynamic variable and NextBinaryState (if any) for which a transition occurs, ie it binary value is 1*)
(* sol1[x_ ]:=FindMin[ListSolsBinary,EqReal,BinaryVars,RealVars,RealStateVars,BinaryStateVars,x,AltBinaryState,NextBinaryState, G[BinaryVars,RealVars,StateVars,r]+Apply[Jstar,NextState]]; *)
(* In sol1[x,y], x is the full state and y is the corresponding solution list from which FindMin2 seeks the one with least cost *)
   sol1[x_,y_List ]:=FindMin2[y,x,AltBinaryState,NextBinaryState, G[BinaryVars,RealVars,StateVars,r]+Apply[Jstar,NextState]];
(* OptSol pairs up states with corresponding solutions and applies sol1 to each pair *)
   OptSol=MapThread[sol1[#1,#2]&,{S,ListSolsAll}];
Print["Finished FindMin. r = "];Print[r];
   Data=Map[Join[S[[#]],{OptSol[[Range[Length[S]],1]][[#]][[1]]}]&,Range[Length[S]]];
   Jstar=Interpolation[Data,InterpolationOrder -> 1];
   ustar=Map[Inner[Rule,BinaryControlVars,Flatten[BinaryControlVars//.OptSol[[#,2]]],List]&,Range[Length[S]]];
(* Print[ustar]; *)
ustarList=Join[ustarList,{ustar}];
];
Simplify[ustarList]
]

Controller Maps`

The following functions are used to generate the lookup table that implements the optimal feedback controller. The main function is BuildLookupTable which opens and closes a textfile in which the lookup arrays are written. BuildLookupTable calls the function LookupTable that actually constructs and writes the arrays.

MyToString[X_List]:=Module[{B,BB},
B=X[[1]];
Map[(B=ToString[B]<>","<>ToString[#])&,Drop[X,1]];
B
]

LookupTable[StateDimensions_,ControVar_,TableName_,Control_,ChannelName_]:=Module[{A,a,B,dRule,A1,Levels,TransLevels,ATest,MatrixNames,channel,ArrayStrings,OnesList},
A=Array[a,StateDimensions];
dRule=MapThread[(#1->ControVar/.#2)&,{Flatten[A],Control}];
A1=A/.dRule;
OnesList=Position[A1,1];
ArrayStrings=Map[(ToString[TableName]<>"("<>MyToString[#]<>") = 1;")&,OnesList];
WriteString[ChannelName,ToString[TableName]<>" = zeros("<>MyToString[StateDimensions]<>");\n"];
Map[WriteString[ChannelName,#<>"\n"]&,ArrayStrings];
]

BuildLookupTable[StateDimensions_List,BinaryControlVars_List,TableNameList_List,Control_,FileName_String]:=Module[{channel},
channel=OpenWrite[FileName];
MapThread[LookupTable[StateDimensions,#1,#2,Control,channel]&,{BinaryControlVars,TableNameList}];
Close[channel]
]

The following is an example.

5. Network with Induction Motor and UPS

Classical Power System

n=number of generators
m=number of PV buses
p=number of PQ buses

Pe[i_,N_]:=(∑_(j=1)^N)(V[i] V[j] G[i,j] Cos[δ[i]-δ[j]]+V[i] V[j] B[i,j] Sin[δ[i]-δ[j]])
Qe[i_,N_]:=(∑_(j=1)^N)(V[i] V[j] G[i,j] Sin[δ[i]-δ[j]]-V[i] V[j] B[i,j] Cos[δ[i]-δ[j]])
f[n_,m_,p_]:=Table[-Pe[i,n+m+p]+P[i],{i,n+m+p}]
g[n_,m_,p_]:=Table[-Qe[i,n+m+p]+Q[i],{i,n+m+1,n+m+p}]
x[n_,m_,p_]:=Join[Table[δ[i],{i,n+m+p}],Table[V[i],{i,n+m+1,n+m+p}]]
LoadFlow[n_,m_,p_]:=Join[f[n,m,p],g[n,m,p]]

Symbolic Model

Clear[B,G];Clear[BB,GG];
BB=Array[B,{2,2}];
GG=Array[G,{2,2}];
For[i=1,i<3,i++,Map[{B[i,#]=0,G[i,#]=0}&,{1,2}]];
B[1,1]=b11=-a/n;B[2,2]=b22=-a/n-d;
B[1,2]=b12=a/n;B[2,1]=b12;G[1,2]=g12=0;G[2,1]=-g12;
G[2,2]=g22=c;
(GG+I  BB)//MatrixForm

FF=f[1,0,1]
GG=g[1,0,1]

Trans={δ[1]->θ[1],δ[2]->θ[2]+θ[1]};

F0=Simplify[FF/.Trans]

G0=Simplify[GG/.Trans]

If V[2] is known

Eqns=Inner[Equal,{F0[[2]],G0[[1]]}/.{Cos[θ[2]] ->A,Sin[θ[2]]->B,P[2]->-P,Q[2]->0},{0,0},List]

Network

Network Solutions

For simplicity take n=1

LogicalExpand[FullSimplify[Reduce[{-P-a B V[1] V[2]-c ((V[2])^2)==0,a A V[1] V[2]-(a+d) ((V[2])^2)==0,A^2+B^2==1,a>0,c>0,d>0,V[1]>0,V[2]>0,P>=0},{A,B,V[1],V[2]},Reals],Reals]]

Notice that there are two possible solution sets. The first corresponds to P=0, the second to P>0.  When P=0, the solution is precisely what we had before. For P>0 the solution structure is more complicated. The next paragraphs consider the two cases: 1) V[1]=2 (saturated field) and 2) V[2]=1 (controlled load bus).

Saturated Field, V[1]=2, P=0

FullSimplify[Solve[B+((A c)/(a+d))==0&&(Sqrt[((((a+d))^2)/((c^2)+(((a+d))^2)))])==A&&P==0&&V[2]==((a A 2)/(a+d)),{A,B,V[2]}],a>0&&c>0&&d>0&&a>d]

Saturated Field, V[1]=2, P>0

SolsA=LogicalExpand[FullSimplify[Reduce[{V[2]==((a A 2)/(a+d)),2==(Sqrt[-(((((a+d))^2) P)/((a^2) A (A c+B (a+d))))]),(Sqrt[1-(A^2)])+B==0,a>0,c>0,d>0,P>0,a>d},{A,B,V[2]},Reals],Reals]]

Length[SolsA]

FullSimplify[ToRadicals[Solve[SolsA[[1]][[{1,2,3}]],{A,B,V[2]}]]]

FullSimplify[ToRadicals[Solve[SolsA[[2]][[{1,2,3}]],{A,B,V[2]}]],a>0&&c>0&&d>0&&P>0]

V2sol=FullSimplify[ToRadicals[Solve[SolsA[[3]][[{1,2,3}]],{A,B,V[2]}]],a>0&&c>0&&d>0&&P>0]

Resolve[Exists[P,4 (a^4)-4 (a^2) c P-(((a+d))^2) (P^2)≥0],Reals]

Reduce[4 (a^4)-4 (a^2) c P-(((a+d))^2) (P^2)≥0&&a>0&&c>0&&d>0,P,Reals]

Notice that 3 solutions were identified. One gives a solution for the bifurcation point. the other 2 produce identical solutions for V[2], one is negative so it is ruled out.

Controlled Load Bus, V[2]=1, P=0

FullSimplify[Solve[B+((A c)/(a+d))==0&&(Sqrt[((((a+d))^2)/((c^2)+(((a+d))^2)))])==A&&P==0&&1==((a A V[1])/(a+d)),{A,B,V[1]}],a>0&&c>0&&d>0]

Controlled Load Bus, V[2]=1, P>0

Suppose V=1

Eqns1=FullSimplify[Reduce[{-P-a B V[1]-c ==0,a A V[1] -(a+d) ==0,A^2+B^2==1,a>0,c>0,d>0,V[1]>0,P>0},{A,B,V[1]},Reals],Reals]

V1sol=FullSimplify[Solve[Eqns1[[{5,6,7}]],{A,B,V[1]}],a>0&&c>0&&d>0&&P>0]

Data

Data={Pm->0.5,Rr->0.25,Rl->2,Xr->0.125,Im->1,ω->2,a->0.375,n->1,Rv->40,Ib->0.04,Cb->1/2}

Set up equations to calculate network bus voltages as required.

c0=((1/Rl) + Rr s/(Rr^2+s^2*Xr^2))/.Data;
d0=((Xr s^2)/(Rr^2+s^2*Xr^2))/.Data;
Data1=Join[Data,{c->(1-η) c0,d->(1-η) d0}]

Network Replacements

Controlled load bus z1 Rule, Saturated field z2 rule.

z1Rule=z1->Chop[Simplify[(((Subscript[δ, qq1])+(Subscript[δ, qq2])+(Subscript[δ, qq3])) ((((Sqrt[(((a+d))^2)+(((c+P))^2)])/a))/.{P->1/Rv})+(Subscript[δ, qq4]) ((((Sqrt[(((a+d))^2)+(((c+P))^2)])/a))/.{P->(1/Rv+Ib)}) + ((Subscript[δ, qq5])+(Subscript[δ, qq6])+(Subscript[δ, qq7]))((Sqrt[(c^2)+(((a+d))^2)])/a))/.Data1]]
z2Rule=z2->Chop[Simplify[(((Subscript[δ, qq1])+(Subscript[δ, qq2])+(Subscript[δ, qq3])) (((Sqrt[((2 (a^2)-c P+(Sqrt[4 (a^4)-4 (a^2) c P-(((a+d))^2) (P^2)]))/((c^2)+(((a+d))^2)))]))/.{P->1/Rv})+(Subscript[δ, qq4]) (((Sqrt[((2 (a^2)-c P+(Sqrt[4 (a^4)-4 (a^2) c P-(((a+d))^2) (P^2)]))/((c^2)+(((a+d))^2)))]))/.{P->(1/Rv+Ib)}) + ((Subscript[δ, qq5])+(Subscript[δ, qq6])+(Subscript[δ, qq7])) ((2 a)/(Sqrt[(c^2)+(((a+d))^2)])))/.Data1]]

const1=((V1>0)∧((-V1)>-2)∧(V2==1))∨(V1==2);
VEqn=((V2==1)⇒(V1==z1))∧((V1==2)⇒(V2==z2));
GenIP[VEqn∧const1,{0≤V2≤2,0≤V1≤2,0≤z1≤2,0≤z2≤2}]

Induction Motor Equivalent Circuit

General Case

(*Parameters:Rs: Stator resistance
Rr: Rotor resistance
RL: Equivalent load resistance
Xm: Magnetizing reactance
Xs: Stator reactance
Xr:  Rotor reactance
Lm: Magnetizing inductance
Vs: Stator voltage
Is: Stator current
I1: Magnetizing current
I2:Rotor current
Is: Stator current 
s: Slip :the relative difference of speeds of field and rotor 
*)

sols=Simplify[Solve[{Is==I1+I2, Vs==(Rs+I Xs+I Xm) I1, Vs==(Rs+Rr/s+I Xs+I Xr) I2},{Is,I1,I2}]]

YY=Is/.sols[[1]]

Y=ComplexExpand[Coefficient[YY,Vs]]

Simplify[Together[(-(Xm/((Rs^2)+(((Xm+Xs))^2)))-(Xs/((Rs^2)+(((Xm+Xs))^2)))-(((s^2) Xr)/((((Rr+Rs s))^2)+(s^2) (((Xr+Xs))^2)))-(((s^2) Xs)/((((Rr+Rs s))^2)+(s^2) (((Xr+Xs))^2))))]]

Simple Case (Rs=0,Xs=0,Xm=∞)

sols=Simplify[Solve[{Is==I2, Vs==(Rr/s+I Xr) I2},{Is,I2}]]

YY=Is/.sols[[1]]

Y=ComplexExpand[Coefficient[YY,Vs]]

s[t]=(ω0-ωm[t])/ω0

D[s[t],t]

Solve[ss==(ω0-ωm)/ω0,ωm]

ωmm[t]=-(-1+ss[t]) ω0;

D[ωmm[t],t]

Vs (((Rr s)/((Rr^2)+(s^2) (Xr^2)))) ((1-s)/s)

Vs Y[[1]] ((1-s)/s)

Simplify[1+(1-s)/s]

Truncate[Vs (((Rr s)/((Rr^2)+(s^2) (Xr^2)))) ((1-s)),s,1]

Truncate[((Rr (1-s) s Vs^2)/((Rr^2)+(s^2) (Xr^2))),{s},1]

?Truncate

V2->1,E->1,c->1,n->1

a->(Sqrt[1+d^2])
a^2->1+d^2

(*   *)
Rr^2+(s^2) Xr^2 ==Rr s

Reduce[Rr^2+(s^2) Xr^2 ==Rr s&&Rr>0&&Xr>0&&0<s<1,{Rr,Xr,s}]

(* Max value for Xr*)
N[(Sqrt[Rr-(Rr^2)])/.Rr->1/4]

(* Value of s *)
N[((Rr/(2 (Xr^2)))-(1/2) (Sqrt[(((Rr^2)-4 (Rr^2) (Xr^2))/(Xr^4))]))/.{Rr->1/4,Xr->0.125}]

(* Value of reactance *)
Rr s/(Rr^2+(s^2) Xr^2 )/.{Rr->1/4,Xr->0.125,s->0.254033}

(* Same value as above  *)
Re[Y]/.{Rr->1/4,Xr->0.125,s->0.254033}

(* Value of suceptance *)
Xr s^2/(Rr^2+(s^2) Xr^2 )/.{Rr->1/4,Xr->0.125,s->0.254033}

(* Same value as above  *)
Abs[Im[Y]]/.{Rr->1/4,Xr->0.125,s->0.254033}

(* Magnitude of the admittance *)
Sqrt[(Rr s/(Rr^2+(s^2) Xr^2 ))^2+(Xr s^2/(Rr^2+(s^2) Xr^2 ))^2]/.{Rr->1/4,Xr->0.125,s->0.254033}

(* Same value as above  *)
Abs[Y]/.{Rr->1/4,Xr->0.125,s->0.254033}

(* Magnitude of the augmented admittance *)
Sqrt[(1/RL+Rr s/(Rr^2+(s^2) Xr^2 ))^2+(Xr s^2/(Rr^2+(s^2) Xr^2 ))^2]/.{Rr->1/4,Xr->0.125,s->0.254033,RL->2}

(* Same value as above  *)
Abs[Y+1/RL]/.{Rr->1/4,Xr->0.125,s->0.254033,RL->2}

Discretization  of the Induction Motor

Pd[s_]:=V2^2 Rr s /(Rr^2+s^2*Xr^2) 
Pe[s_]:=Pd[s] (1-s)

(* solution=DSolve[s'[x]-((Pm-V1^2 s[x](1-  Rr s[x]))/(Rr^2+s[x]^2 Xr^2))/(Im ω^2)==0,s[x],x] *)

Solve[((Pm-V2^2*Rr*s*(1-s)/(Rr^2+s^2*Xr^2))/.Join[Data1,{V2->1}])==0,s]

F[x_]:=((x*0.9442-V2^2*Rr*s*(1-s)/(Rr^2+s^2*Xr^2))/.Join[Data1,{V2->1}]);
SS=s/.Transpose[Map[Solve[F[#]==0,s]&,Range[944]/944]];

SS1=MapThread[{#1,#2}&,{Range[944]/944,Flatten[SS[[1]]]}];
SS2=MapThread[{#1,#2}&,{Range[944]/944,Flatten[SS[[2]]]}];

gr3=ListPlot[SS1,PlotJoined->True];
gr4=ListPlot[SS2,PlotJoined->True,PlotStyle->GrayLevel[0.75]];
Show[gr3,gr4,AxesLabel->{"(Subscript[P, m])","s"}];

gr1=ListPlot[Flatten[SS[[1]]],PlotJoined->True];gr2=ListPlot[Flatten[SS[[2]]],PlotJoined->True,PlotStyle->GrayLevel[0.75]];
Show[gr1,gr2,AxesLabel->{"(Subscript[P, m])","s"}];

Options[ListPlot]

f[x_]:=(1-η)(Pm-Pe[x])/(Im ω^2);
g[x_]:=FullSimplify[-x+x0+(f[x0]+f[x])h/2]
J[x_]:=Evaluate[FullSimplify[D[g[x],x]]]

g[x]

FullSimplify[J[x]]

FullSimplify[g[x0]]
FullSimplify[J[x0]]

SlipEq=Chop[(FullSimplify[Together[x0-g[x0]/J[x0]]]/.x0->s)]/.Data/.h->0.5

Battery

ChargeEq=σ+((1/Cb) (h (((Subscript[δ, qq1])+(Subscript[δ, qq2])+(Subscript[δ, qq3])+(Subscript[δ, qq7])) 0 + (Subscript[δ, qq4]) Ib - ((Subscript[δ, qq5])+(Subscript[δ, qq6])) /Rv)))/.Data

6. Optimal Control

Logical Specification

B3=exactly[1,{q1,q2,q3,q4,q5,q6,q7}]∧exactly[1,{qq1,qq2,qq3,qq4,qq5,qq6,qq7}]∧
((q1∧(¬s1))⇒(qq2))∧((q1∧c)⇒qq4)∧((q1∧(¬((¬s1)∨(c))))⇒(qq1))∧
((q2∧(¬s2))⇒(qq3))∧((q2∧s1)⇒(qq1))∧((q2∧d)⇒qq5)∧((q2∧(¬(s1∨(¬s2)∨d)))⇒qq2)∧
((q3∧s2)⇒qq2)∧((q3∧d)⇒qq6)∧((q3∧¬(s2∨d))⇒qq3)∧
((q4∧(¬c))⇒qq1)∧((q4∧(¬((¬c))))⇒qq4)∧
((q5∧(¬d))⇒qq2)∧((q5∧(¬s2))⇒qq6)∧((q5∧(¬((¬d)∨(¬s2))))⇒qq5)∧
((q6∧s2)⇒qq5)∧((q6∧(¬d))⇒qq3)∧((q6∧(¬(s2∨(¬d))))⇒qq6)∧
(q7⇒qq7)

AAA=GenIP[B3,{s1,s2,c,d,q1,q2,q3,q4,q5,q6,q7,qq1,qq2,qq3,qq4,qq5,qq6,qq7,0≤σ≤1}]

GenIP[q1⇒¬s2∧¬d∧atmost[1,{c,¬s1}],{q1,s1,s2,c,d}]

GenIP[q1⇒(¬s2∧¬d∧atmost[1,{c,¬s1}]),{q1,s1,s2,c,d}]

ControlLogic={
(q1⇒(¬s2∧¬d∧atmost[1,{c,¬s1}]∧¬qq3∧¬qq5∧¬qq6∧¬qq7∧exactly[1,{qq1,qq2,qq4}])),
(q2⇒(¬c∧atmost[1,{d,¬s2,s1}]∧¬qq4∧¬qq6∧¬qq7∧exactly[1,{qq1,qq2,qq3,qq5}])),
(q3⇒(¬s1∧¬c∧atmost[1,{d,s2}]∧¬qq1∧¬qq4∧¬qq5∧¬qq7∧exactly[1,{qq2,qq3,qq6}])),
(q4⇒(¬s1∧¬s2∧¬d∧¬qq2∧¬qq3∧¬qq5∧¬qq6∧¬qq7∧exactly[1,{qq1,qq4}])),
(q5⇒(¬s1∧¬c∧atmost[1,{¬d,¬s2}]∧¬qq1∧¬qq3∧¬qq7∧exactly[1,{qq2,qq5,qq6}])),
(q6⇒(¬s1∧¬c∧atmost[1,{¬d,s2}]∧¬qq2∧¬qq1∧¬qq4∧exactly[1,{qq5,qq6,qq3,qq7}])),
(q7⇒(¬c∧¬d∧¬s1∧¬s2)∧¬qq1∧¬qq2∧¬qq3∧¬qq4∧¬qq5∧¬qq6∧exactly[1,{qq7}])};

BBB=GenIP[ControlLogic[[1]],{s1,s2,c,d,q1,q2,q3,q4,q5,q6,q7,qq1,qq2,qq3,qq4,qq5,qq6,qq7,0≤σ≤1}]

StateLogic={(q1⇒¬(qq3∨qq5∨qq6∨qq7)),
(q2⇒¬(qq4∨qq6∨qq7)),
(q3⇒¬(qq1∨qq4∨qq5∨qq7)),
(q4⇒¬(qq2∨qq3∨qq5∨qq6∨qq7)),
(q5⇒¬(qq1∨qq4∨qq3∨qq7)),
(q6⇒¬(qq2∨qq1∨qq4))};

BinaryStateVars={(Subscript[δ, q1]),(Subscript[δ, q2]),(Subscript[δ, q3]),(Subscript[δ, q4]),(Subscript[δ, q5]),(Subscript[δ, q6]),(Subscript[δ, q7])};
Delta=Array[del,Length[BinaryStateVars]];
Map[(del[#]=0)&,Range[Length[Delta]]];
StateReplaceList=Map[Inner[Rule,BinaryStateVars,ReplacePart[Delta,1,#],List]&,Range[Length[BinaryStateVars]]]

ControlReplacements=ControlRules[ControlLogic,StateReplaceList,{q1,q2,q3,q4,q5,q6,q7,qq1,qq2,qq3,qq4,qq5,qq6,qq7},{c,d,s1,s2}];

ControlReplacements[[1]]

ControlReplacements[[2]]

ControlReplacements[[3]]

ControlReplacements[[4]]

ControlReplacements[[5]]

ControlReplacements[[6]]

ControlReplacements[[7]]

Map[Dimensions[#]&,ControlReplacements]

Setup Functions

const1=((V1>0)∧((-V1)>-2)∧(V2==1))∨(V1==2);
A7=B3∧VEqn∧const1;
AA=GenIP[A7,{s1,s2,c,d,q1,q2,q3,q4,q5,q6,q7,qq1,qq2,qq3,qq4,qq5,qq6,qq7,0≤σ≤1,0≤η≤1,0≤V2≤2,0≤V1≤2,0≤s≤1,0≤z1≤15,0≤z2≤10}]

eqBinary=BinaryEquations[AA,{V1,V2,s,σ,z1,z2,η}]

varBinary=BinaryVariables[AA,{V1,V2,s,σ,z1,z2,η}]

Length[varBinary]

eq1=(eqBinary/.{(Subscript[δ, q1])->1,(Subscript[δ, q2])->0,(Subscript[δ, q3])->0,(Subscript[δ, q4])->0,(Subscript[δ, q5])->0,(Subscript[δ, q6])->0,(Subscript[δ, q7])->0})/.ControlReplacements[[1,5]]
Position[eq1,False]

DDD=Map[Position[(eqBinary/.{(Subscript[δ, q1])->1,(Subscript[δ, q2])->0,(Subscript[δ, q3])->0,(Subscript[δ, q4])->0,(Subscript[δ, q5])->0,(Subscript[δ, q6])->0,(Subscript[δ, q7])->0})/.ControlReplacements[[1,#]],False]&,Range[Length[ControlReplacements[[1]]]]]

FFF=Position[DDD,{}]

ControlReplacements[[1]][[Flatten[FFF]]]

LongList=Timing[Reduce[eq1,Union[Flatten[Map[Variables[Level[#,1]]&,eq1]]],Integers]];

LongList[[1]]

Length[LongList[[2]]]

eqReal0=Complement[AA,eqBinary]

eqReal=(eqReal0/.{z1Rule,z2Rule})/.{η->((Subscript[δ, qq2])+(Subscript[δ, qq5])) 0.4 + ((Subscript[δ, qq3]) + (Subscript[δ, qq6]) + (Subscript[δ, qq7])) 0.8};

(eqReal)/.{(Subscript[δ, qq1])->0,(Subscript[δ, qq2])->0,(Subscript[δ, qq3])->0,(Subscript[δ, qq4])->1,(Subscript[δ, qq5])->0,(Subscript[δ, qq6])->0,(Subscript[δ, qq7])->0,s->0.5}

eqReal0

((z1/.z1Rule)/.{η->((Subscript[δ, qq2])+(Subscript[δ, qq5])) 0.4 + ((Subscript[δ, qq3]) + (Subscript[δ, qq6]) + (Subscript[δ, qq7])) 0.8})/.{(Subscript[δ, qq1])->1,(Subscript[δ, qq2])->0,(Subscript[δ, qq3])->0,(Subscript[δ, qq4])->0,(Subscript[δ, qq5])->0,(Subscript[δ, qq6])->0,(Subscript[δ, qq7])->0,s->0.5}

((z2/.z2Rule)/.{η->((Subscript[δ, qq2])+(Subscript[δ, qq5])) 0.4 + ((Subscript[δ, qq3]) + (Subscript[δ, qq6]) + (Subscript[δ, qq7])) 0.8})/.{(Subscript[δ, qq1])->1,(Subscript[δ, qq2])->0,(Subscript[δ, qq3])->0,(Subscript[δ, qq4])->0,(Subscript[δ, qq5])->0,(Subscript[δ, qq6])->0,(Subscript[δ, qq7])->0,s->0.5}

BinTest=BinarySolutions[eqBinary,StateReplaceList,ControlReplacements];

Length[BinTest[[3]]]

Length[BinTest[[3]][[1]]]

TransitionPossibilities=Union[Map[Drop[#,9]&,BinTest[[3]]]]

Union[Map[Drop[#,9]&,BinTest[[1]]]]

Union[Map[Drop[#,9]&,BinTest[[2]]]]

{{(Subscript[δ, c])->0,(Subscript[δ, d])->0,(Subscript[δ, qq1])->0,(Subscript[δ, qq2])->0,(Subscript[δ, qq3])->1,(Subscript[δ, qq4])->0,(Subscript[δ, qq5])->0,(Subscript[δ, qq6])->0,(Subscript[δ, qq7])->0,(Subscript[δ, s1])->0,(Subscript[δ, s2])->0},{(Subscript[δ, c])->0,(Subscript[δ, d])->0,(Subscript[δ, qq1])->0,(Subscript[δ, qq2])->1,(Subscript[δ, qq3])->0,(Subscript[δ, qq4])->0,(Subscript[δ, qq5])->0,(Subscript[δ, qq6])->0,(Subscript[δ, qq7])->0,(Subscript[δ, s1])->0,(Subscript[δ, s2])->1},{(Subscript[δ, c])->0,(Subscript[δ, d])->0,(Subscript[δ, qq1])->1,(Subscript[δ, qq2])->0,(Subscript[δ, qq3])->0,(Subscript[δ, qq4])->0,(Subscript[δ, qq5])->0,(Subscript[δ, qq6])->0,(Subscript[δ, qq7])->0,(Subscript[δ, s1])->1,(Subscript[δ, s2])->1},{(Subscript[δ, c])->0,(Subscript[δ, d])->1,(Subscript[δ, qq1])->0,(Subscript[δ, qq2])->0,(Subscript[δ, qq3])->0,(Subscript[δ, qq4])->0,(Subscript[δ, qq5])->1,(Subscript[δ, qq6])->0,(Subscript[δ, qq7])->0,(Subscript[δ, s1])->0,(Subscript[δ, s2])->1}}

Union[Map[Drop[#,9]&,BinTest[[4]]]]

Union[Map[Drop[#,9]&,BinTest[[5]]]]

{{(Subscript[δ, c])->0,(Subscript[δ, d])->0,(Subscript[δ, qq1])->0,(Subscript[δ, qq2])->1,(Subscript[δ, qq3])->0,(Subscript[δ, qq4])->0,(Subscript[δ, qq5])->0,(Subscript[δ, qq6])->0,(Subscript[δ, qq7])->0,(Subscript[δ, s1])->0,(Subscript[δ, s2])->1},{(Subscript[δ, c])->0,(Subscript[δ, d])->1,(Subscript[δ, qq1])->0,(Subscript[δ, qq2])->0,(Subscript[δ, qq3])->0,(Subscript[δ, qq4])->0,(Subscript[δ, qq5])->0,(Subscript[δ, qq6])->1,(Subscript[δ, qq7])->0,(Subscript[δ, s1])->0,(Subscript[δ, s2])->0},{(Subscript[δ, c])->0,(Subscript[δ, d])->1,(Subscript[δ, qq1])->0,(Subscript[δ, qq2])->0,(Subscript[δ, qq3])->0,(Subscript[δ, qq4])->0,(Subscript[δ, qq5])->1,(Subscript[δ, qq6])->0,(Subscript[δ, qq7])->0,(Subscript[δ, s1])->0,(Subscript[δ, s2])->1}}

Union[Map[Drop[#,9]&,BinTest[[6]]]]

{{(Subscript[δ, c])->0,(Subscript[δ, d])->0,(Subscript[δ, qq1])->0,(Subscript[δ, qq2])->0,(Subscript[δ, qq3])->1,(Subscript[δ, qq4])->0,(Subscript[δ, qq5])->0,(Subscript[δ, qq6])->0,(Subscript[δ, qq7])->0,(Subscript[δ, s1])->0,(Subscript[δ, s2])->0},{(Subscript[δ, c])->0,(Subscript[δ, d])->1,(Subscript[δ, qq1])->0,(Subscript[δ, qq2])->0,(Subscript[δ, qq3])->0,(Subscript[δ, qq4])->0,(Subscript[δ, qq5])->0,(Subscript[δ, qq6])->1,(Subscript[δ, qq7])->0,(Subscript[δ, s1])->0,(Subscript[δ, s2])->0},{(Subscript[δ, c])->0,(Subscript[δ, d])->1,(Subscript[δ, qq1])->0,(Subscript[δ, qq2])->0,(Subscript[δ, qq3])->0,(Subscript[δ, qq4])->0,(Subscript[δ, qq5])->1,(Subscript[δ, qq6])->0,(Subscript[δ, qq7])->0,(Subscript[δ, s1])->0,(Subscript[δ, s2])->1}}

Union[Map[Drop[#,9]&,BinTest[[7]]]]

Example

Data1

ContinuousS=Flatten[Outer[List,Range[0.1,0.5,0.1],Range[0.25,1,0.25]],1];
h=0.5;
m=50;
DiscreteS=Range[7];
S=Map[Flatten[#]&,Flatten[Outer[List,ContinuousS,DiscreteS,1],1]];
r=1/25;
G[BinaryVars_List,RealVars_List,StateVars_List,n_Integer]:=(((V2-1)^2+(σ-1)^2+r η^2+10  (Subscript[δ, qq7]))/m)/.{η->((Subscript[δ, qq2])+(Subscript[δ, qq5])) 0.4 + ((Subscript[δ, qq3]) + (Subscript[δ, qq6]) + (Subscript[δ, qq7])) 0.8};
GN[StateVars_List]:=0;  (* terminal cost *)
BinaryStateVars={(Subscript[δ, q1]),(Subscript[δ, q2]),(Subscript[δ, q3]),(Subscript[δ, q4]),(Subscript[δ, q5]),(Subscript[δ, q6]),(Subscript[δ, q7])};
BooleanStateVars={q1,q2,q3,q4,q5,q6,q7,qq1,qq2,qq3,qq4,qq5,qq6,qq7};
DiscreteState=BooleanStateVars;
NextBinaryState={(Subscript[δ, qq1]),(Subscript[δ, qq2]),(Subscript[δ, qq3]),(Subscript[δ, qq4]),(Subscript[δ, qq5]),(Subscript[δ, qq6]),(Subscript[δ, qq7])};
BinaryControlVars={(Subscript[δ, c]),(Subscript[δ, d]),(Subscript[δ, s1]),(Subscript[δ, s2])};
BooleanControlVars={c,d,s1,s2};
RealStateVars={s,σ};
RealVars={V1,V2};
AllRealVars={V1,V2,s,σ};
BinaryVars=Complement[varBinary,BinaryStateVars];
DynamicsRHS={SlipEq,ChargeEq}/.{η->((Subscript[δ, qq2])+(Subscript[δ, qq5])) 0.4 + ((Subscript[δ, qq3]) + (Subscript[δ, qq6]) + (Subscript[δ, qq7])) 0.8};
{T,Control}=Timing[OptimalPolicy2[eqBinary,eqReal,ControlLogic,BinaryVars,RealVars,BinaryStateVars,BooleanStateVars,NextBinaryState,RealStateVars,DynamicsRHS,BinaryControlVars,BooleanControlVars,S,GN,G,m]];

Control[[-1]]

Control[[-1]]-Control[[-2]]

T

S

Length[S]

Lookup Table

SetDirectory["C:\Projects\PowerSystems\STTR Power System Management\Notebooks\Power System Examples\UPS"]

BuildLookupTable[{5,4,7},BinaryControlVars,{"AAA","BBB","CCC","DDD"},Control[[-1]],"TableDefs.m"]

8. References

[1] Li, Q., Y. Guo, et al. (1999). Transformation of Logical Specification into IP-formulas. 3rd International Mathematica Symposium (IMS '99), Hagenburg, Austria, Computational Mechanics Publications, WIT Press.
[2] Li, Q., Y. Guo, et al. (2000). "Modelling Integer Programming with Logic: Language and Implementation." IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences E83-A (8): 1673-1680.
[3] Hooker, J. (2000). Logic-Based Methods for Optimization: Combining Optimization and Constraint Satisfaction, Wiley-Interscience.
[4] Bemporad, A. and M. Morari (1999). "Control of Systems Integrating Logic, Dynamics, and Constraints." Automatica 35: 402-427.
[4] Tyler, M. L. and M. Morari (1999). "Propositional Logic in Control and Monitoring Problems." Automatica 35 (4): 565-582.
[5] Williams, H. P. (1993). Model Building in Mathematical Programming, John Wiley and Sons.
[6] Matthew G. Earl, Raffaello D'Andrea, Multi-vehicle Cooperative Control Using Integer Programming
[7] J. E. Beasley ." Advances in Linear and Integer programming", Oxford Science Publication, 1996
[8] H. G. Kwatny, E. Mensah, D. Niebur, G. Bajpai, C. Teolis, “Logic Based Design of Optimal Reconfiguration Strategies for Ship Power Systems,” Proceedings 7th IFAC Symposium on Nonlinear Control Systems, August 22-24, 2007, Pretoria, South Africa.